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SUMMARY 


A  new  family  of  Fourier-based  estimators  of  the  parameters 
of  the  multivariate  Gaussian  distribution  is  presented.  The  estimators 
are  equivalent  to  parametric  density  estimators.  Three  distinct 
estimators  arise,  each  of  which  is  robust  and  reduces  to  the  maximum 
likelihood  estimator  as  a  special  case.  By  varying  the  window  width 
of  a  parametric  density  estimator,  a  set  of  diagnostics  which  are 
useful  in  problems  of  outlier  detection  and  clustering  are  obtained. 

An  example,  using  a  trivariate  data  set,  is  given.  >  . 
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1.  INTRODUCTION 


In  fundamental  papers  Rosenblatt  (1956),  Parzen  (1962),  and 
Watson  and  Leadbetter  (1963)  considered  the  problem  of  estimating  a 
density.  While  the  approach  used  by  these  authors  is  non-parametric , 
we  shall  be  concerned  with  density  estimation  in  a  parametric  frame¬ 
work.  In  this  framework  the  form  of  the  density  is  assumed  to  be  known 
apart  from  parameters  which  completely  define  the  density.  We  shall 
exclusively  consider  the  problem  of  estimating  the  parameters  of  the 
p-variate  normal  density 

f(x)  =  , 2"D j  iexp(-4(x  -  u)TD  1(x  -  u))  ,  (1.1) 


where  U  is  the  mean  vector  and  D  is  the  covariance  or  dispersion 
matrix. 

The  characteristic  function  corresponding  to  (1.1)  is  given 
by 


5(u)  =  exp(i-.Tu  -  ‘j  uTDu)  (1.2) 

X 

where  u  =  (u. ,  u« ,  . . . ,  u  )  is  a  column  vector  of  real  numbers.  If 

~  i  /  p 

x^,  x^,  ...»  :<n  is  a  random  sample  from  (1.1),  the  sample  characteris¬ 
tic  function 

-1  -  T 

0  (u)  *  n  )  exp  ( iu  x.)  (1.3) 

n  -  ^  -  -J 


is  unbiased  for  (1.2)  for  every  fixed  u  .  Information  concerning 

and  D  can  be  extracted  from  :  (u)  bv  methods  similar  to  those  of 
~  n  - 

Paulson  and  Nicklin  (193  3).  In  this  case  the  estimators  of  u  and  D 
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would  be  determined  from  score  functions  for  u  and  D  developed  by 
different iating 


15:' 


|$(u)  “  in(u) j  2  exp(-m“u^Ku)du 


(1.4) 


with  respect  to  y  and  D  ,  setting  K  ■  D  and  then  setting  the  re¬ 
sulting  expressions  to  zero.  The  expression  (1.4)  is  therefore  not  a 
bona  fide  objective  function. 

An  alternate  approach  is  to  proceed  as  follows.  Let 


)  (u)  ■  n  1  V  exp(iuTD  ^(x .  -  - 

a  jii  *  -  ^  - 


>)  . 


(1.5) 


Whereas  $n(u)  is  a  statistic,  ;n(u)  is  not,  but 


E(-;n(u))  -  V (u)  -  exp(-l5u  u) 


(1.6) 


The  parameters  -  and  D  of  (1.5)  may  be  estimated  by  making  '^(u) 
and  ^(u)  match  up  in  some  sense.  There  are  many  ways  in  which  this 
may  be  done,  but  we  shall  present  only  that  which  has  been  found  to  be 
most  theoretically  as  well  as  computationally  convenient  and  useful. 
Computational  convenience  is  not  a  matter  of  small  importance,  it  is 
in  fact  the  overriding  consideration  when  the  dimension  p  becomes 
large;  an  estimator  is  effectively  worthless  if  it  cannot  be  calculated. 

We  shall  ultimately  obtain  three  pairs  of  estimators  if 
(u,  D)  ,  all  arising  from  consideration  of  matching  up  .^(u)  with 
'l(u)  .  All  pairs  of  estimators  will  be  useful  la  the  identification  of 
potential  outliers,  in  exploratory  data  analysis,  and  in  clustering 
problems.  Furthermore,  ail  estimators  will  te  qualitatively  robust  in 


WWW 


the  sense  that  their  multidimensional  influence  functions  will  be 
bounded  and  re-descendant;  all  will' reduce  to  familiar  estimators 
a  special  case. 
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2.  DERIVATION  OF  THE  ESTIMATORS  BASED  ON  "J  (u) 

n  ~ 


Define  Che  function 


Qm  n(y*  D> 

m,  n  ~ 


9  9  T 

i*(u)  -  'P  (u)  |  exp(-m“u  u)du  .  (2.1) 


The  function  Q  (U,  D)  is  a  bona  fide  obiective  function  for  ’  and 
m,n  -  —  -  - 

D  when  0  <  tn.  < 00  .  The  integrand  of  equation  (2.1)  can  be  rewritten 
in  terms  of  the  residuals  in  u  , 


R(u)  *  exp(-!i(l  +  m“)uTu)  -  n  ^  £  exp(iu^D  ^(x .  -  u)  -  4m^u^u)  (2.2) 


j-1 


-J 


as 


Qm  n(U’  D) 
m,  n  ~ 


R(u)  R  (u)du  * 


lR(u) 1  du  , 


where  R  (u)  denotes  the  complex  conjugate  of  R(u)  . 

Estimators  which  result  from  minimizing  (2.1)  are  thus  those 
which  minimize  the  integrated  squared  moduli  of  residuals  defined  by  the 
Fourier  transforms  in  (2.2). 

Explicit  integration  of  (2.1)  yields 

_;5Pr  -2.  2.-4p  -  -  ■  2.-1  T 

[n  (m  )  ■  exp\ -(-n  )  (v  -  v,  )  (v.  -  v,  )  ; 

r  r  -  j  '-x  - ; 

J  x 

-  2  n  ^(m^  +  :.)  /  exD'-(2  +  im~)  ‘v‘  v  ■ 

7  -J  j 


+  (m“  +  1 ) 


'  -Pi 


i  2.3) 


t 


Equating  the  partial  derivatives  of  (2.3)  with  respect  to  y  and  D 
(Dwyer,  1967)  to  zero  gives  the  estimating  equations  for  y  and  D  as 
respectively 


n  ^7r'‘P(m^  +  ^)  ^  ^  D  ^(x.  -  a  )exp'.-(4m“  +  2)  ^(x.  -  u)TD  2(x.  -  y): 

j  -J  -  ~  '  ~J 

(2.4) 


and 


TJSP[.?Sn-2(m2)”,sP"1  l  l  D"L(x  -  x  )(x  -  x  )TD 

j  k  ~J  -k  ' 


T.-l 


exp  {-(4  m2)  L(x^  -  x^1  D  1(x^  -  xk) } 


15n'1(m2  +  's)”’iP"1  7  D_1(x.  -  y 


l  “  -  :)(?j  - '^Y1 

j  J 


;xp  {-  (4  m2  +  2)_1  (x .  -  u)T  D~L(x.  -„);]=  0  . 


(2.3) 


From  (2.4)  and  (2.5)  we  find  that  the  estimators  for  „  and  D  , 


and  D 


i  » 


say,  satisfy  the  implicit  equations 
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where 


A  -  (m2  +  ‘-5)  ^P  1  i  (x.  -  _)(:<.  -  -)T 
J  -J 


x  exp-. -(Am“  +  2)  L(x.  -  - )^D  L(x.  -  ~):  , 


(2.8) 


B  »  (2n)  2  (m")  ‘P  2  >  >  (x.  -  x.  )(x.  -  x,  ) 

-r  r  -J  -«i  ~j  -x 
J  « 


exp\-(Am“)  2(x .  -  x.  )  TD  2  (x  .  -  x,  ) ;  , 

-j  ~j  ~<c 


(2.9) 


vljm  "  exP  HA  m2  +  2)  i(Xj  -  2)T  D  L(Xj  -  u) } 


(2.10) 


T  -1 

The  generalized  distances  (x.  -  -)  D  (x.  -  . )  and 
T  -1 

(x^  -  x^)  D  (x^  -  x^)  constitute  an  essential  component  of  the  esti¬ 
mation  process  defined  by  (2 .6)-(2 .10)  .  These  distances  have  been  used 
in  a  variety  of  contexts  including  clustering  (Gnanades ikan ,  1977,  Ch . 

A),  identification  of  potential  outliers  (Barnett  and  Lewis,  1978,  Ch . 

6;  Cnanades ikan  and  Kettenring,  1972;  Rohlf,  1973),  and  regression 
diagnostics  (Belsley,  Kuh,  and  Welsch,  1980,  Ch .  2).  Because  the  diag¬ 
nostic  capability  of  these  generalized  distances  is  embodied  in  the 
estimation  process  itself,  we  have  in  a  real  sense  combined  the  processes 
of  estimation  and  criticism.  The  word  criticism  is  used  here  in  the 
sense  of  Daniel  (1959),  Box  (1979,  p.  2),  and  Paulson  and  Nicklin  (1983). 
Criticism  is  the  process  of  assessing  the  internal  '  'ns  iscen -v  _■>:  tne 
data  and  the  tentative  model,  in  this  case  that  the  x.  are  indepen¬ 
dently  distributed  according  to  a  p-variate  Gaussian  law. 


As  m  -►  00  ,  the  estimator  -*■  n  A  xj  >  the  ordinary 
sample  mean  vector.  The  covariance  matrix  D  is  not  strictly  estimable 

A 

through  (2.7)  when  m  00  since  we  end  up  with,  at  p  =  y  ,  the 
identity 


i  (;j  -  -  “i> 


~~  i,  i  l  %  -  JkXjj  -  ;k> 


(2.11) 


which  does  not  explicitly  involve  D  .  Either  side  of  (2.11)  may  be 
used  as  an  estimator  for  D  on  the  truth  of  the  independent,  identical, 

/V  A 

p-variate  normality  of  the  x^  .  The  estimators  and  are  not 

m-estimators  because  of  their  explicit  dependence  on  the  pairwise 
generalized  distances.  The  comment  that  either  side  of  (2.11)  may  be 
used  to  produce  estimators  for  D  leads  us  to  presently  develop  two 
different  estimators  for  D  . 


3.  SIMPLIFIED  ESTIMATORS  OF  D 


The  estimator  is  numerically  difficult  to  extract  from 

(2.7) ,  especially  if  p  is  large.  Fixed  point  and  Newton-Raphson 

algorithms  for  producing  have  been  developed  and  found  to  work  well 

except  for  computational  expense.  However,  the  intriguing  identity 

(2.7)  at  the  solution  ^  and  suggest  two  alternative  estimators 

for  D  .  The  first  is  based  on  A  of  (2.8).  The  joint  estimators  for 
U  and  D  are  obtained  by  determining  the  value  in 

S  *  7  {K  (x.  -  u)(x.  -  U)T  -  D;exp'-(4n^  +  2)  ^(x.  -  _)TD  ^Cx.  -  _); 

~2  j  2  ~j  *•  ~J  -  -J  ~  J  ~ 

(3.1) 

which  makes  £(3^)  =  0  .  Direct  computation  gives 


2  +  2m 


1  +  2m 


(3.2) 


Thus  the  joint  estimators  for  u  and  D  ,  and  6,  ,  satisfy  the 

implicit  equations 


=  1 


!  X  .  V,  . 

h  J 


L  v^r 


(3.3) 


T* 

0  l  (x.  -  U)(X  -  „)  "v 

2  +  2m"  i 


1  +  2m" 


L 


J]m 


(3.4) 


with 


exp: -(4m“  *  1< 


..) " d~" r x .  -  .  )' 


(3.5) 


Another  pair  of  estimators  for  u  and  D  may  be  based  on 
B  of  the  right  hand  side  of  equation  (2.8).  We  require  the  constant 


in 


S3  -  l  l  CK3(Xj  -  xk) (Xj  -  xk)L  -  Dl  exp{-4m2)  1(x.  -  xk>ip~i(xj  -  xR)' 
j  k 
j^k 

(3.6) 


.  T  -1 , 


which  makes  EtS^)  =  0  .  Since  x.  -  x^  ,  j  j*  k  ,  is  11^(0,  2D)  we 


~3  ~  -l  ~k 

find  by  direct  computation  that 


K 


3  -  TorrTT  «  ♦  -?>  +  V5  ■ 


(3.7) 


m 


Thus  the  joint  estimators  for  y  and  D  ,  y^  and  63  ,  satisfy  the 
implicit  equations 


y  y  x  vv. 

U  -j  3jkn 


n 


(3.8) 


j*k 


3jkm 


”  ~k)(^j  -  ^k5  V  3  jkn 


D  =  K. 


Jiii 


(3.9) 


jA 


3  jkm 


with 


1  - L  T  -1 

v  =  exp-fln")  (x,  -x.)  D  (x.  -  x,  ) 
3jkm  -j  ~k  ~j  -k 


(3.10) 


Both  estimators  D?  and  63  are  computationally  more 
attractive  than  63  .  All  of  6  ,  j  =  1,  2,  3  ,  are  positive  definite 
with  probability  one  for  m  >  0  whenever  n  >  p  although  all  may  be 
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Figure  1  provides  some  of  these  efficiencies  as  a  function 

A  /A 

of  m  and  p  •  The  efficiencies  associated  with  (y^  has  the 

A  .  A  ■ 

same  asymptotic  efficiency  as  )  and  D3  are  approximately 

/S  A 

those  of  and  °2  •  A  small  Monte  Carlo  study  suggests  that  the 

efficiencies  related  to  are  slightly  higher  than  those  of  • 

The  efficiencies  of  the  off-diagonal  element  of  relative 

to  the  maximum  likelihood  estimates  of  covariances  are  very  nearly  equal 
to  those  of  the  diagonal  elements. 

The  contours  of  all  influence  functions  of  the  estimators 

A  A 

y j  ,  Dj  at  the  multivariate  normal  are  all  closed  and  bounded.  This 
closedness  property  implies  that  estimation  may  be  combined  with 
clustering  or  may  be  used  as  a  clustering  algorithm  or  used  to  evaluate 
the  results  of  a  clustering  algorithm.  The  clustering  capability  asso¬ 
ciated  with  estimation  of  y  and  D  allows  for  an  identification  of 
potential  outliers  in  multivariate  normal  data. 
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4.  PARAMETRIC  DENSITY  ESTIMATION 

A  A 

The  estimators  and  of  Section  2  are  equivalent  to 

those  derived  from  the  following  parametric  density  estimation  consid¬ 
erations.  The  expression  R(u)  of  (2.2)  represents  a  difference  of 
characteristic  functions  whose  inverse  is 


r(x)  -  g(x)  -  gn(x)  (4.1) 

where  g(x)  is  the  spherical  normal  distribution  with 

g(x)  »  i  2rr  ( 1  +  m^)lj  ^exp.-Hd  +  m“)  \x~x:  (4.2) 

and  gn(x)  is  an  estimator  of  g(x)  with 

g  (x)  ■  n_1  7  '2~m“I  ^exp- -(.2m  )  lzTz.-  ,  (4.3) 

5n  -/  ’  ~  ~J  •  } 

Zj  ■  (x  -  D  5(x.  -  '-))  .  (4.4) 


The  matrix  I  is  the  pxp  identity  matrix.  The  expression  g^(x)  is 
unbiased  for  g(x)  when  the  x^  are  p-variate  Gaussian.  By  the  multi¬ 
dimensional  version  of  Parseval's  theorem 

r  ?  r  o 

Q_  (y,  D)  ■  ' R(u) ' *du  ■  (2-)?  R"(x)dx  (4.5) 

'n,n  Jr  ~  '  •  R 

P  P 

Accordingly,  the  estimators  and  are  those  arrived  at  from 

minimizing  the  integrated  squared  residual  in  x  and  may  be  accurately 
termed  parametric  density  estimators. 

The  robustness  characteristics  of  the  estimators  and 

are  easily  seen  from  working  with  (4.5),  e.g.,  differentiation  of  (4.5) 
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with  respect  to  y  gives 
f  3  g  (x) 

— ;■  (g(x)  -  g  (x))dx  -  0  . 

rp 

From  this  equation  we  readily  deduce  for  0  <  m  <  00  that  the  marginal 

contribution  of  x^  ,  say,  to  the  estimate  of  y  depends  on  all  the 

Xj  ,  not  just  xr  ,  as  well  as  the  assumed  p-variate  Gaussian  density, 

that  it  is  bounded,  and  that  it  is  redescendent  as  llx  !!  ■+■  <®  .  The 

n 

function  gR(x)  is  similar  to  a  Parzen  kernel  density  estimator;  the 

2  T 

form  of  the  kernel  is  determined  by  the  function  exp(-m  u  u)  of 
(2.1).  The  advantage  to  the  choice  of  the  Gaussian  kernel  is  that  it 
leads  to  convenient  computational  expressions.  Choice  of  functions 
[oj(u)j  other  than  exp(-rrTu  u)  generally  leads  to  kernels  which 
do  not  permit  closed  form  integration  of  (4.5)  and  thus  make  for  in¬ 
tractable  numerics  when  p  >  2  .  Intractable  numerics  result  from  the 
need  to  compute  integrals  by  numerical  methods. 

Indeed,  the  potential  for  application  of  parametric  density 
estimation  procedures  is  bleak  if  integrals  of  the  form  (2.1)  or  (4.5) 
for  a  non-Gaussian  p-variate  density  f(x)  cannot  be  evaluated  in 
closed  form.  The  reason  is  simple:  high  order  numerical  integration 
for  p  large  is  not  viable  even  for  today's  computers. 
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5.  AN  EXAMPLE  AND  DISCUSSION 


A  choice  of  m  is  needed  to  implement  the  estimators  of 
Sections  2  and  3.  We  feel  that  exploration  of  the  data  vls-a-vis  the 
p-variate  normal  model  is  in  many  cases  more  useful  since  it  leads  to 
a  deeper  understanding  of  the  data  generating  process  and  the  actual 
problem  at  hand.  Accordingly,  it  will  often  be  useful  to  use  a  range 
of  values  of  m  in  our  analyses  since  the  responses  of  the  parameter 
estimates  and  final  weights  to  changes  in  m  generates  very  useful 
diagnostics . 

Table  1  gives  the  chemical  analysis  of  20  geological  speci¬ 
mens  in  terms  of  percentage  of  iron  (Fe),  sodium  (Na),  and  potassium 
(K)  which  had  been  believed  by  geologists  to  be  approximately  homo¬ 
geneous  in  character  and  also  trivariate  Gaussian.  All  two  way  scatter- 
plots  of  this  data  are  given  in  Figure  2.  A  few  points  visually  "stick 
out"  enough  to  question  the  belief  of  homogeneity  and  simple  3-variate 
Gaussianity.  It  is  of  interest,  however,  to  determine  how  the  three 
robust  estimators  of  covariance  react  to  this  data  for  various  values 
of  m  . 

Table  2  presents  the  estimates  of  the  variances  and  correla- 

2  2 
tions  for  m  *  10,  5,  2  and  maximum  likelihood.  At  m  ■  10  ,  all 

estimates  are  reasonably  close  to  those  of  maximum  likelihood  and  at 

2  2 
m  *5  there  are  seme  minor  differences.  However,  at  m"  =  2  we  find 

that  there  are  now  important  differences  in  the  estimates.  These  dif¬ 
ferences  provide  a  warning  that  something  may  be  awry  regarding  the 
internal  consistency  of  the  data  and  the  assumption  of  3-variate  normal¬ 
ity.  Of  particular  interest  are  the  estimates  of  the  correlation 
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TABLE  1 

i 

CHEMISTRY  OF  GEOLOGICAL  SPECIMENS 

M 

Observation 

Number  Fe  Na 

y/i 

A  2.6  1.7 

1 

34  ^ 

B  2.1  2.1 

3.4 

C  1.3  2.8 

I-7  '  P 

D  2.2  2.1 

3.0  V 

v% 

E  1.3  2.2 

's' 

3-6  & 

F  1.6  2.2 

3.7  i 

G  3.1  2.2 

2.4 

H  2.8  1.7 

3.8  ■•*>' 

I  4.0  1.4 

3-3  S 

J  2.6  1.9 

.  ?  ‘"'i* 

3.2  X 

K  1.5  2.0 

4-2  i 

L  3.9  1.6 

25  i 

M  3.1  1.7 

3  5 

N  3.1  1.9 

3-6  $ 

0  1.7  2.1 

3.6  ■« 

sg 

P  1.4  2.2 

4-0  $$ 

Q  3.0  2.0 

3.8  | 

R  2.9  2.2 

5.5  fli 

S  2.9  2.8 

v!y 

5.9 

•\  * 

T  2.9  2.9 

6 . 5  ;-.;j 

i 

1 

V*  1 
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FIGURE  2b.  Percentages  of  iron  and  potassium  for  geological  specimens 
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for  the  four  estimators.  Maximum  likelihood  perceives  a  relatively 
strong  positive  correlation  between  Na  and  K  ,  the  estimators  D, 
and  perceives  a  relatively  weak  positive  correlation  between  Na 

and  K  ,  and  che  estimator  perceives  a  relatively  strong  negative 

correlation  between  Na  and  K  .  An  examination  of  Figure  lc  leads 
one  co  conclude  that  each  interpretation  is  permissible,  depending  on 
which  data  points  are  considered  inconsistent  with  the  remainder  of 
the  data  and  the  3-variate  Gaussianity.  It  should  be  emphasized  that 
a  purely  robust  procedure  which  focused  on  95’i  efficiency  may  not 
have  picked  up  a  potential  problem  with  the  data-model  unit  while  an 
exploratory  (sensitivity  of  solutions  to  changes  in  m  )  analysis  did. 
Of  course,  for  this  3-dimensional  situation  we  can  easily  see  by 
graphical  methods  that  the  data  are  not  homogeneous  and  that  there  may 
be  some  interesting  geological  structure  which  warrants  further  inves¬ 
tigation.  Since  simple  graphical  methods  may  prove  less  effective  in 
higher  dimensions,  there  is  a  real  advantage  to  using  these  estimators 
in  combination  with  graphical  methods  for  exploration  and  deeper  under¬ 
standing  of  data  and  its  generation. 


The  final  weights  v .  .  and  w  . .  ,  i.e., 

°  •  i  m  ’<m 


eva.ua tea 


at  ,  D,  ,  are  useful  in  diagnostic  examination  of  the  data.  Table 

3  presents  the  final  observational  weights  for  m “  *  2  ,  v  for 

J  ••* 

1  *  1,  2  and  /  v  for  D .  .  All  three  estimation  procedures 

k*j  3jKr'  '3 

call  attention  to  observations  18,  19,  and  20  as  being  potentially 
different  from  the  remainder  in  the  assumption  of  Guassianity.  Methods 

2  and  3  strongly  indicate  that  observation  3  is  potentially  different 
while  the  estimator  1  weakly  indicates  che  potential  difference,  ie  do 
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TABLE  3 


FINAL  OBSERVATIONAL  WEIGHTS  FOR  ra' -  2 


Dhflprvfl  t  ”f  nn 

H 
<  Q 

°2 

A 

°3 

Number 

2  , 
m  =  2 

m2  -  2 

m2  -  2 

1 

0.058 

0.061 

0.059 

2 

0.069 

0.078 

0.074 

3 

0.038 

0.005 

0.011 

4 

0.068 

0.069 

0.066 

5 

0.062 

0.067 

0.064 

6 

0.065 

0.073 

0.070 

7 

0.051 

0.033 

0.034 

8 

0.067 

0.069 

0.066 

9 

0.054 

0.052 

0.051 

10 

0.068 

0.077 

0.074 

11 

0.060 

0.062 

0.059 

12 

0.049 

0.048 

0.047 

13 

0.067 

0.073 

0.070 

14 

0.056 

0.055 

0.055 

15 

0.067 

0.075 

0.072 

16 

0.060 

0.067 

0.064 

17 

0.042 

0.036 

0.038 

18 

0.6*10'3 

0.5*10_4 

0.007 

19 

0.1*10-6 

0.3*10"9 

0.010 

20 

l.lxlO-8 

-12 

0. 3*10 

0.010 
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